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Abstract 



A recently developed quasi two-dimensional exact-exchange formalism within 
the framework of Density Functional Theory has been applied to a strongly 
inhomogeneous interacting electron gas, and the results were compared with 
state-of-the-art Variational Quantum Monte Carlo (VMC) numerical simula- 
tions for a three-dimensional electron gas under a strong external potential. 
The VMC results, extremely demanding from the computational point of 
view, could be considered as a benchmark for the present theory. We observe 
a remarkable qualitative and quantitative agreement between both methods 
from the comparison of the exchange- hole densities, exchange-energy densi- 
ties, and total exchange-energies per particle. This agreement is increasingly 
improved with the strength of the external potential when the electron gas 
becomes quasi-two-dimensional. 
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The search for accurate calculation schemes within the framework of Density Functional 
Theory (DFT) is of paramount importance in atomic, molecular, and solid state physics, [l 
While DFT maps the problem of interacting inhomogeneous electron systems to a problem 
of effectively non-interacting inhomogeneous electrons, it provides no help in the determi- 
nation of the exchange-correlation component of the total energy functional. |2|,|3| Largely, 
this obstacle has been overcome by using the so-called Local Density Approximation (LDA), 
which is based in the approximate equivalence in the long- wavelength limit between the inho- 
mogeneous and homogeneous interacting electron gases, with the latter taken as a reference 
system. || The LDA (and generalized gradient corrections to the LDA, as the GGA) has 
been widely employed, however, even for strongly inhomogeneous systems, such as atoms 
and molecules, with surprisingly good results. || The application to strongly inhomoge- 
neous solid state systems has been equally successful, || but some warnings raised quite 
recently when applied to two-dimensional electron gases (2DEG), @ of the type formed at 
the interface between two semiconductors such as GaAs/Al x Gai_ x As. || Recently a study 
by Nekovee, Foulkes, and Needs (NFN) using Variational Quantum Monte Carlo (VMC) 
numerical simulations showed significant failures of the LDA and GGA in the description 
of the exchange hole of a 3D electron gas under the effect of a strong periodic potential. || 
Since plane wave expansions are a routine approach in DFT ab-initio codes NFN results call 
for methods that could describe accurately the electron gas in regions where the chemical 
bond is formed. 

In this paper we present the first application of a recently developed quasi two- 
dimensional exact-exchange DFT (XX-DFT) method to a strongly inhomogeneous "metal- 
lic" system; we have found that this local exact exchange theory overcomes the problems 
presented by LDA and GGA. Previous application of the XX-DFT formalism includes atomic 
and molecular systems, [[TtJ semiconductors (insulators), |TT| and quasi 2DEG's. The 
method is distinguished from LDA, by avoiding the approximation of taking the exchange 
energy of the inhomogeneous electron gas as given by the exchange energy of the homoge- 
neous electron gas, evaluated at the same (local) value of the density of the inhomogeneous 
system. Instead, the XX-DFT formalism uses the fact that the exchange energy is an ex- 
plicit function of the occupied orbitals, |L3[ and exploit this to derive a local potential that 



minimizes a Kohn-Sham energy which includes the exchange energy exactly. |L4| This elab- 
orate treatment of the exchange potential, beyond the standard LDA-DFT scheme, was 
shown to be crucial in providing a theoretical explanation of striking experimental results in 
asymmetric semiconductor quantum wells. [|TjJ These experiments and theoretical analysis 
suggest that the occupancy of a quantum well subbands proceed in an abrupt way, mediated 
by inter-subband exchange, the explanation being quite natural in the XX-DFT scheme. 

The system under study is shown schematically in Fig. 1, and has been motivated by 
the recent analysis of a similar system by NFN using VMC numerical simulations; f|] the 



same system has been also studied by Rushton, Tozer, and Clark, |15| using the so-called 
Weighted Density Approximation (WDA), that treats both exchange and correlation in an 
approximate, but non-local, way. |TB| Our model consists of a cosenoidal double well external 
potential V ext (z) = Vocos(g,z), which strongly modulates the electronic density along the z 
direction, but keeps the translational symmetry in the x — y plane; the system is ultimately 
confined by two infinite barriers located at z = ± {d + A). The charge neutrality of the 
system was fulfilled by means of a positive jellium slab of uniform density n = 



s ■ 
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where r s = ro/ao, with Atttq/S being the volume per electron, and ao = 0.529 A the Bohr 
radius. Following Ref. (9), r s was fixed at r s = 2, while the amplitude Vo was chosen to 
be Vo = 2.08e° F , where e F is the Fermi energy corresponding to n; these values corresponds 
approximately to those of Al. We studied systems with three different modulations q of 
the external potential, corresponding to q\ = l.llA;^, q 2 = 1.55k F , q 3 = 2.17 k F , with 
o-okp = (S^n) 1 / 3 . Note that in our model, decreasing q leads to an increasing A, and 
consequently to increasingly smooth potentials (as the potential barrier height Vo is fixed), 
approaching the three-dimensional uniform situation. 

Before starting with calculations, it is important to realize that A = 2ir/q ~ 1 / k F ; this 
leads to a strong quantization of the electronic degrees of freedom along z. Following the 
terminology of the quasi-2DEG's, we will denote as "subbands" each one of these quantized 
discrete levels. As electrons are free to move in the x — y plane, an in-plane parabolic 
dispersion relation is associated to each one of these subbands. To proceed, and exploiting 
the translational symmetry along the x — y plane (area A), we propose as solutions of 
the three-dimensional Kohn-Sham (KS) equations (j)^ (r) = exp(i k • p)£° (z) /yA, with 
k the in-plane wave vector, p the in-plane coordinate, (z) the subband wave function 
corresponding to an electron with a spin projection a (j or j) and subband index v {y — 
0, 1, 2, ...). Proceeding with these solutions, the KS equations take the form of effective one- 
dimensional Schrodinger-like equations [[!]]] 
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2dz 2 



+ V KS {z, a) 



£{z) = el£{z), (1) 



with Vks( z , a) being the local KS spin-dependent potential, and the eigenvalues. Vks(z, a) 
is the sum of several contributions 

V KS (z, a) = V ext {z) + V H (z) + V£(z, a) + V c (z, a), (2) 

with Vjj, V*, and V c being the Hartree, exchange, and correlation contributions, respectively. 
The index % stands for the exact-exchange {% = XX) or Local Density Approximation 
(i = LDA). The Hartree potential is given by its classical expression, 

X+d 

V H {z) = -2-k J dz' \z - z'\ [n(z') - n jeU (z)} , (3) 

-(A+d) 

with nj e u(z) = n 9(\ + z)9(X — z), and 9(x) being the step function (9(x) — 1 if x > 0, 
9{x) = if x < ). For V c (z, a) we follow the DFT prescription and define it through the 
functional derivative with respect to the density of the correlation contribution E c to the 
total energy, 

V c (z,a) = -^- (4) 
dn(z, a) 

and use LDA for its practical evaluation, using a recent parametrization of the correlation 
energy per particle for the interacting homogeneous electron gas by Gori-Giorgi, Sacchetti, 
and Bachelet. |18| In equations above, n(z) is the zero-temperature three-dimensional elec- 
tron density 
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n ( z )= Yl 



n[z, a 



7? 



(5) 



n(z, a) being the fraction of a polarized electrons. \x is the chemical potential (or Fermi level) 
of the system, which is determined by the electrostatic and thermodynamic equilibrium with 
a reservoir, and k v p = y2{ji — e%)9(/i— e%). It remains to define V xx (z). As explained above, 
we avoid its usual local (LDA), or semilocal (GGA) expressions, using instead the recently 
developed XX-DFT formalism. TO, Tl] As shown elsewhere, [|I^] the exact-exchange potential 
for our quasi two-dimensional geometry is defined by the equation 



V x XX (z,a) 
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SV KS (z',a) 
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(6) 



5n (z, a) 



In equation above, the functional derivatives SE X /S^ (z) and SE x /5k F a can be evaluated 
directly from the explicit expression for the exchange energy E x in terms of (z) and 
k v p . (z) /5Vks( z ' \ a ) an d Sk F a /5Vks ( Z ,< J ) are evaluated using first-order perturbation 
theory from Eq. (1). Finally, x~ x [z, z') = SVks ( z> , °") l& n is the inverse of the operator 
Xa(z, z') = 5n (z,a) /5V KS {z', a) given by 



Xtr(z, z') 
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The first term in Eq. (|7]) comes from first-order perturbation theory [|Tl|], whereas the 
second term results from first-order perturbation theory and the thermodynamic equilibrium 
between the electronic system and the reservoir that fixes a common chemical potential 
H allowing (in principle) the change of the number of particles. Indeed, without the 
reservoir, the operator Xo-{z,z') in general cannot be inverted because it is singular ||TT 
The alternative to V* x (z) is the widely employed V^ DA (z), defined as 



V x LDA (z,a) = - 



Qn(z, a) 



71 



1/3 



(8) 



As it is well known, Eq.(l), together with Eqs.(3), (4), and (6) for XX-DFT, or (7) for LDA- 
DFT should be solved iteratively in a self-consistent way. All results denoted as XX (LDA) 
in what follows have been obtained by performing this numerical self-consistent procedure 
with V* x (z) {V x DA {z)) inserted in Eq.(l), and by assuming a paramagnetic situation for the 
present high-density metallic system. Using as input the set of parameters defined above, the 
output of this self-consistent calculation correspond to a situation with two (four, counting 
spin) occupied subbands {y = 0, 1) for the three different modulations q/k F . 

One way to characterize the accuracy of a many-body approach, is to analyze the 
exchange- correlation hole density, in terms of which the exchange-correlation energy could 
be defined. As the emphasis of the present work is on exchange effects, we will restrict our 
analysis mainly to the exchange-hole density. We start from the exact expression for the 
exchange-hole density, [[2(1] 
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Mr , r+R)s _iw^, (9) 

where 

p(r, r + R) = J2 D* - <( k )l (* + R) </w (r) (10) 

is the density matrix associated to our problem. e£(k) = e a v + (ftk) 2 /2m and note that 
p(r, r) = n(r). The physical meaning of h x (r, r + R) is that it represents the density of the 
exchange hole at point r + R (observational point) due to the presence of an electron (test 
particle) located at r. Evaluating Eq.@ in our quasi two-dimensional geometry of Fig. 1, 
we obtain 

h x {r, r + R) = -—-J—— KkiJxik^Jxik^CAz + Z)t v {z%>(z + Z)£,(z), 

(11) 



where R = (R n , Z), and Ji(x) is the first-order cylindrical Bessel function. |21j Taking the 



limit R n , Z — > 0, it is easy to check that h x (r, r) = — n(r)/2, and by explicit integration the 



fulfillment of the exact constraint [BO 



dR/ix(r,r + R) = -1, (12) 

is immediate. In Eq.(|ll]), it is interesting to note, first, the absence of the in-plane coordinate 
p corresponding to the test particle , which result from the translational symmetry in the 
x — y plane; choosing p = without any loss in generality, h x (r, r + R) could be equally 
written as h x (z;z + Z,R n ). Second, the dependence of h x only on the magnitude of R n 
means that the exchange hole has cylindrical symmetry along the z axis. Third, integration 
of n(z) x h x (z; z + Z,R n )/R (R = |R|) in z,Z, and R u leads exactly to the definition of 
the exchange energy used to derive V x xx (z) in Eq.(6). In the limit kpR u ^> 1 (assuming 
a finite occupancy of all occupied subbands), and using the asymptotic expansion of the 
Bessel function, we obtain 

h x (r,r + R) : 



(Tri^Mr) 

y/w4 cos(£^R M - ^) cos^'i?,, - ^)Q(z + Z)Uz)i v ,{z + Z)&(z). (13) 



From Eq. ([T3|) , is evident that the asymptotic behavior of h x (r, r + R) is given by R~ 3 along 
the x—y plane; on the other side, the decay along z is controlled by the extent of the subband 
wave-functions £ v (z). In what follows, we will denote by h xx (r,r + R) the exchange-hole 
density as given by Eq.([ll]), evaluated with the self-consistent subband wavefunctions £„(z) 
that are the solutions of Eq.(l) with V xx (z). On the other side, evaluation of Eq.@ in the 
LDA-DFT scheme using 3D plane waves yields 
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h™\v,v + R) = -?f- F[k F (v)R], (14) 

with F{x) = 9 [(sinx — xcosx)/x 3 ] 2 . In Eg . flli]) , and according to the LDA prescription 
kp(r) = [37r 2 n(r)] 1 ^ 3 . We will denote by h^ DA (r, r + R) the exchange-hole density as given 
by Eq.(|TJ|), evaluated with the self-consistent subband wavefunctions £ v (z) that are the 
solutions of Eq.(l) with V£ DA {z). It is easy to see that F(x) — > 1 for x — > 0, while it decays 
as x~ 4 for x = kF(r)R 3> 1. This simple dependence of F(x) on x implies that h^ DA is 
always centered on the electron coordinate r. Besides, it is clear that h% DA (r,r + R) is a 
spherically symmetric function of R, and centered at r. 

Fig. 2 allows a qualitative comparison of both holes at a density maximum, corresponding 
to the electron coordinate p = 0, z = —A/2 being located at the left potential well of Fig. 1. 
For this case, both holes are centered at the electron coordinate . Note that for both upper 
and lower panels, only half of the exchange hole has been represented; the full hole could 
be obtained by reflection with respect to the z axis. The exact-exchange hole is contracted 
along the confinement direction, due to the fact allowed above that the decay along the x — y 
plane and along the z direction are different, as seen explicitly from Eq. (|l"3"D. The fact that 
h xx (r,r + R) has essentially zero strength at the right well minima, simply means that 
the system is in the "atomic" or isolated well limit, with a very small amount of tunneling 
or subband wave-function overlap between the two wells. It is also interesting to note the 
Friedel-like oscillations of the exact-exchange and LDA holes, each of them reflecting the 
different types of symmetry of the corresponding holes. 

The situation displayed in Fig. 3, corresponding to the electron coordinate p = 0, z = 
at the density minimum, is strikingly different, compared with Fig. 2. Here, h x (r,r + 
R) attains its maximum strength when the observational coordinate Z coincides with the 
density maximum (z = ±A/2). This is the physically correct behavior that one expects 
for the exchange- hole density, and that can easily be obtained from Eq. (|TTj) : Restricting for 
simplicity the analysis to the case where only one subband [y = 0) is occupied, this equation 
simplifies to 

h xx (r, r + R) = - lJl(fc ;f" )]2 Uz + Z)f , (15) 

which clearly attains its maximum strength, as a function of Z, when the argument of the 
subband wave function [z + Z) coincides with a density maximum. This non-local behavior 
of the exact-exchange hole is clearly preserved when more than one subband is occupied 
(two subbands are occupied for the situation displayed in Figs. 2, 3, and 4). For z = ±A/2 
(Fig. 2) this happens for Z = 0, while for z = 0, this happens for Z = ±A/2. The LDA 
exchange hole, on the other side, is centered at the electron coordinate z. Attention should 
be paid also to the quantitative values of h xx (r, r + R) and h% DA (r, r + R) in Fig. 3. In 
the exact-exchange case, h x has a maximum strength value of about half of its value in 
Fig. 2, which is reasonable keeping in mind the constraint given by Eq. (|12[). On the other 
side, h^ DA has a maximum strength of about two orders of magnitude smaller than h xx , 
and accordingly it has an "anomalous" large extent in real space in order of satisfy Eq. (|12|). 
Note also that the extent in real space for h xx (r,r + R) is similar in Figs. 2 and 3, as 
it should be, as in this case and for both situations the size of the hole is dictated by the 
accumulation points of electronic density. 
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We show in Fig. 4 the exchange-energy densities differences, defined as follows: 
Ae^ MC (z) = e xx {z)-e v x MC {z) (full thick line), and Ae L x DA (z) = e x x (z) - e L x DA {z) (dashed 
thick line). Exchange-energy densities are obtained from 

400 = 5 J ^n(r)/4(r,r + R) (16) 

with i = XX, LDA. The full thin line corresponds to the XX density n xx (z), the dotted 
thin line to the VMC density n VMC (z). We will discuss first the comparison between XX 
and VMC calculations. While in principle both methods give an exact treatment to the 
exchange-energy, and accordingly they should coincide, in practice there are three possible 
sources for the small differences observed in Fig. 4: i) while our double quantum well system 
is a very good approximation to the periodic system studied by NFN, the equivalence is not 
perfect from the "model" point of view, ii) even if both models were identical, for the 
practical implementation of the VMC calculation, the external potential V ex t(z) used by NFN 
is not identical of ours (although it should be close to it), and Hi) even if both the model 
and the external potential were identical, still we are treating correlation at the LDA level, 
while NFN treats correlation with VMC accuracy, presumably a much better approximation 
than LDA, and definitively quite different. As a consequence of these three factors, small 
differences between the two different densities profiles result, which in turn translate to the 
small differences for the exchange-energy densities observed in Fig. 4. It is important to 
realize, however, that points i and ii above are just consequence of the different scheme of 
calculations, while point Hi is the important and intrinsic difference between XX-DFT and 
VMC. It is also interesting to note that the differences between XX and VMC are mainly 
concentrated in regions of high density, as expected; besides, differences increase from top 
(small q) to bottom (large q), as both models increasingly disagree in the large q limit. 
While from the results presented in Table I for the corresponding integrated magnitudes 
the difference between XX and VMC is smallest for the intermediate case q/k° F = 1.55, 
from Fig. 4 it is clear that this is due to an (accidental) error cancellation. Based in the 
almost perfect agreement between n xx (z) and n VMC (z) displayed in Fig. 4 for q/k F = 1.11, 
it is clear that our modelization of the periodic system studied by NFN is at best in this 
quasi two-dimensional limit (see below). It is then plausible to assume that the remaining 
discrepancy between XX and VMC is mainly due to the different treatment of correlation 
in both cases. In this context, it is interesting to point out that in our case, the correlation 
energy per particle is typically 15 % of the exchange-correlation energy per particle. The 
comparison between XX and LDA exchange-energy densities shows large local deviations 
between the two methods, the amplitude of the differences being about a factor of two 
greater than between XX and VMC. 

The situation is similar for the differences in the integrated energies presented in Table 
I, the differences between XX and VMC being much smaller than between XX and LDA for 
q/k F = 1.11 and 1.57. The difference between XX and LDA exchange energies (per particle) 
is smaller than between XX and VMC for q/k F = 2.17. As can be seen from Fig. 4, this 
is in part due to a large error cancellation of the LDA results, which, being non-systematic 
is difficult to control, and also due to the fact that as q increases, our model and the NFN 
model move towards different limits (see below). 

We show in Fig. 5 the exchange energy per particle e x , calculated in different approxi- 
mations, as a function of the potential modulation q/k° F . In order to place these results in 



7 



perspective, it is important to note that the exchange energy per particle for the homoge- 
neous 3D electron gas is 



while for a 2D homogeneous electron gas it is 



^ = (18) 



s 



For r s = 2, we obtain e 3 x D (r s = 2) ~ -0.23, e 2 x D (r s = 2) ~ -0.30. For q/k% -> 0, the 
system resembles a three-dimensional uniform electron gas, which explains the decreasing 
(in absolute value) of the exchange energy for q/k° F < 1. While our XX-DFT formalism is 
ideally suited for the study of strongly inhomogeneous systems, it is numerically ineficcient 
for reaching the 3D limit, when a large number of subbands become occupied and the 
(numerical) calculation becomes increasingly cumbersome. [2"3| The tendency of the exchange 



energy towards the homogeneous 3D limit e^P(2) ~ —0.23 is however clearly displayed in 
Fig. 5. It is also interesting to note how LDA becomes an increasingly better approximation 
to XX for q/kp — > 0, as one expects. The opposite limit q/k F > 1 is also worth analyzing. 
In this case, the periodic model studied by NFN becomes increasingly three-dimensional, 
as the fixed-height barriers becomes thinner and electrons are increasingly able of tunnel 
through them. Our model, on the other side, in the limit q/k F > 1, and as a consequence 
of the constraint A = 2ir/q, departs from the NFN model and becomes more similar to an 
infinite barrier quantum well of width 2d with a positively charged layer at its center. The 
different limits reached by the model of NFN and ours explains the progressive disagreement 
between XX and VMC which is seen for q/k F > 2. Besides these two limiting features, 
clearly the more interesting region is the intermediate regime about q/k F ~ 1, where the 
inhomogeneity reaches its maximum strength. The proximity of the e^ x and e^ MC to the 
2D value e^ D {2) ~ —0.30 tells us that for q/k F ~ 1 the system is quite close to the strict 
two-dimensional limit. Quite good agreement is observed between XX and VMC in this 
strongly inhomogeneous regime, as according to Table I, the percentage error between both 
calculations is about 0.65 %. This should be compared with the corresponding error between 
XX and (our) LDA, that reaches a value of about 2.4 % in the same quasi two-dimensional 
regime, almost four times higher. 

In summary, we have checked the accuracy of the exact-exchange density functional 
method as applied to a strongly inhomogeneous interacting electron gas, and compared 
with the best available results. We have analyzed the exact-exchange hole density, and found 
that it is strongly non-local when evaluated in regions of low-density. From the comparison 
with the results of Variational Quantum Monte Carlo, we found a remaining discrepancy 
of about 0.65 % between the exact-exchange and VMC exchange contributions, that we 
attribute to the different treatment of correlation in both approaches. Exact-Exchange 
Density Functional Theory works at best in the worse case for LDA, that is, when the 
strongly inhomogeneous electron gas is closer to the quasi two-dimensional limit. Our results 
motivate the urge for better and more accurate correlation functionals, and work in this 
direction is in progress. 
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FIGURES 



FIG. 1. Schematic representation of the external potential V ex t(z) = Vocos(qz). The 
wave-length of the modulation is A = 2n/q, and the amplitude Vq = 2.08 e F . d = 2 ao — 1 
A is the distance between the end of the jellium slab (shaded region) and the infinite potential 
barriers located at z = ±(<i + A). 

FIG. 2. Exchange-hole density at a density maximum ( z = —A/2 in Fig. 1), for q/k F = 1.11. 
Upper panel: h xx (r, r + R). Lower panel: h^ DA (r, r + R). 

FIG. 3. Exchange-hole density at a density minimum ( z = in Fig. 1), for q/k° F = 1.11. 
Upper panel: h xx (r,r + R). Lower panel: h^ DA (r, r + R). It is important to realize that the 
numerical value of h% DA at its point of maximum strength coincides with the value of h xx at its 
point of minimum strenght, marked with a square full dot. 



FIG. 4. Left-scale, full (dashed) thick line: Difference between XX and VMC (LDA) ex- 
change-energy densities Ae^ MC (z) = e xx (z) — e^ MC (z) (e xx (z) — e^ DA (z)) , versus position in 
the system, and for the three modulations. Right-scale: exact-exchange (full thin line) and VMC 
(dotted thin line) electronic densities n xx (z) and n VMC (z), respectively; note that both sets of 
densities are almost indistinguishable except for q/k F = 2.17and \z\ — A. Data corresponding to 
el MC {z) and n VMC (z) provided by the authors of Ref.(23). 



FIG. 5. Exchange-energy per particle as a function of the potential modulation. Full and 
dashed lines are our XX (e xx ) and LDA {£x DA ) results, respectively. Discrete points are VMC 
data from Ref. 9. 
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TABLES 



TABLE I. Exchange energy (hartrees) per particle from three different methods: XX-DFT 
(e* x ), LDA-DFT (e% DA ) , and VMC {e^ MC , from Ref. 23). 



q/k° F 


F xx 


P LDA 


P VMC 
^x 


F XX p VMC 
^x ^x 


F XX F LDA 

^x ^x 


1.11 


-0.2911 


-0.2795 


-0.2930 


0.0019 


-0.0116 


1.55 


-0.2739 


-0.2687 


-0.2756 


0.0017 


-0.0052 


2.17 


-0.2469 


-0.2508 


-0.2534 


0.0065 


+0.0039 
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